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Abstract. We study various box-size scaling techniques to obtain the multifractal properties, in terms of 
the singularity spectrum /(a), of the critical eigenstates at the metal-insulator transition within the 3-D 
Anderson model of localisation. The typical and ensemble averaged scaling laws of the generalised inverse 
participation ratios are considered. In pursuit of a numerical optimisation of the box-scaling technique we 
discuss different box-partitioning schemes including cubic and non-cubic boxes, use of periodic boundary 
conditions to enlarge the system and single and multiple origins for the partitioning grid are also imple- 
mented. We show that the numerically most reliable method is to divide a system of linear size L equally 
into cubic boxes of size / for which L/l is an integer. This method is the least numerically expensive while 
having a good reliability. 
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1 Introduction 

The multifractal analysis (MFA) [1-3] of critical electronic 
wavefunction probabilities IV'iP in the three-dimensional 
Anderson model of localisation [4-7] is based on the scal- 
ing of the generalised inverse participation ratio (gIPR) 
Pq. For a system with volume L'^, the gIPR for normalised 
wavefunctions is defined as Pq{l) = X]fe' Mfc(0 where the 

integrated measure //fc(/) = ^^ IV'jj is computed in all Ni 
boxes with linear size I covering the system. At criticality 
the scaling law 

■(.) (1) 



sis. In the present work we analyse the performance of sev- 
eral partitioning schemes, some of which have previously 
been used in the literature to perform the MFA [3,8, 10] . In 
particular, we shall study the application of cubic versus 
non-cubic boxes, the use of periodic boundary conditions 
to enlarge the system, and multiple origins for the parti- 
tioning grid, and adaptive linear fits. In Fig. [T] we indicate 
these various strategies schematically. Here, we show that 
the use of cubic boxes with integer L/l seems to be best 
suited to calculate the singularity spectrum using the box- 
size scaling of Pq . 



Pg(A) CC X^^ 



is expected to hold in a certain range of values for A = l/L 
[1]. The multifractal character is most often studied in 
terms of the singularity spectrum f{a). Its true behaviour 
at criticality can only be found in the thermodynamic 
limit. This limit can be reached in ([T|) by taking the box 
size to be / -^ 0. In implementing this so-called box-size 
scaling method, one usually considers a fixed system size L 
and partitions it into smaller boxes. The scaling behaviour 
of the gIPR with box size is then obtained by varying / and 
averaging over many samples. Hence, with only one sys- 
tem size to be considered, box-size scaling is numerically 
inexpensive and has been much used previously [3,8,9]. 

There exist, however, multiple ways of carrying out 
this box partitioning, and some of them might lead to bet- 
ter results, e.g. better fits or an improved statistical analy- 
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In Refs. [11, 12] we had shown that the box-size scaling 
method is more likely to be affected by finite size effects 
as compared to the so-called system-size scaling approach. 
The latter method is based on varying system size L with 
fixed box size I in the limit L — > oo [13]. Nevertheless, we 
emphasize that box-size scaling is still a useful alternative 
since (i) the computational and data-storage requirements 
are much less demanding than for system-size scaling, (ii) 
in certain regions of the spectrum such as close to the 
maximum, box-size scaling is also quite accurate, and (iii) 
in some situations it might be the only applicable method, 
i.e. when carrying out an MFA of experimental data [14- 
16] for which the system size cannot be easily changed. It 
is therefore important to know how to maximise the use 
and performance of box-size scaling. 
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(b) Unrestricted values of L In 
this case L = 6, Z = 4 and for 
the scaling A — l/L' 
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(c) Non-cubic boxes with 
L = 6, Ix = 2 and ly = 
3. For the scaling A = 
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(d) Multiple origins for cubic boxes and integer L/l. The I non-equivalent origins for the 
partitioning arc shown for the case L = 6 and / = 2 

Fig. 1: Two-dimensional illustration of the different box-partitioning sehemes. The black dashed lines mark the parti- 
tioning grid. The thick grey line corresponds to the boundaries of the physical system with linear size L. In cases (b) 
and (d), the empty points in the shaded regions outside the system are obtained using periodic boundary conditions 
to properly complete all boxes. The position of the black sites highlights the periodicity pattern. 



2 The Anderson model of localisation and its 
numerical diagonalisation 

We use the tight-binding Anderson Hamiltonian in lattice 
site basis as given by 



^-E^ 






XjI, 



(2) 



where site i = {x,y, z) is the position of an electron in a 
cubic lattice of volume V = L^, tij are nearest- neighbour 
hopping amplitudes and e^ is the random site potential 
energy. The hopping amplitude is taken to be i = 1 and 
periodic boundary conditions (PBC) are used to minimise 
boundary effects. We consider Si to have a box probability 
distribution in the interval [— Wc/2, Wc/2], where Wc is 
taken to be the strength of the critical value of the dis- 
order. We assume Wc = 16.5, above which all eigenstates 
are localised [17-20]. We have considered eigenstates ip = 
XizV'il*) only in the vicinity of the band centre E = 0. 
Moreover, we only take about five eigenstates in a small 
energy window at -E = for any given realization of disor- 
der [11, 12]. Since the aim of the present study is to com- 
pare various box-partitioning schemes, we shall restrict 
ourselves to moderate system sizes of L^ = 60'^ = 216000. 
The resulting sparse matrices have been diagonalised us- 
ing J AD AMILU [21]. 



3 Typical and ensemble averaging for the 
MFA 

The well known singularity spectrum f{a) is defined from 
the T{q) exponents via a Legendre transformation: fq = 
f{aq) = qaq — T{q) and aq = T'{q). The physical mean- 
ing of the /(a) is as follows. It is the fractal dimension 
of the set of points where the wavefunction intensity is 
\ipi\'^ ^ L^", that is in a discrete system the number of 
such points N^ scales as L^^"\ The numerical multifractal 
analysis is based on an averaged form of the scaling law ([T]) 
for the gIPR in the limit A = l/L — > 0, where the contri- 
butions from several finite-size critical wavefunctions cor- 
responding to different disorder realizations are properly 
taken into account. In the following, we give the typical- 
and ensemble-averaged form of the singularity strength a 
and singularity spectrum /(a). The expressions are writ- 
ten to be optimal for direct numerical computation. Their 
compact forms are given in Refs. [11, 12]. The scaling law 
for the typical average of the moments Pq is defined as 



e('nP,(A))^;^r'-(g) 



(3) 



where (• • • ) denotes the arithmetic average over many real- 
izations of disorder, i.e. over many different wavefunctions 
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at criticality. The scaling exponents are then defined by 

(lnF,(A)) 



T*yP(q) = lim 



A^o In A 



(4) 



Applying a Legendre transformation [1], we obtain the 
definitions for a'^P and f^P, 

q J q J 

«r=l-^(^|.?.(A)lnM.(A)^ (5a) 



^^<'^-}^.-x 



Nx 



n^T^E/^wi^^^^w 



- anP,(A)) 



(5b) 



The scaling law for the ensemble average involves the 
arithmetic average of Pq over all realizations of disorder, 

(P,(A))(xA^°"°(«). (6) 

Thus the definition of the scaling exponents is 

ln(P,(A)) 



r™''(g) = lim 



A^o In A 



(7) 



and the corresponding definitions of a™*^ and /°"^ can be 



written as 



Nx 



lim ■ -— - 

A^olnA(P,(A)) 



^M?.(A)lnA*fe(A) (8a) 



^K-) = p-oha 



q 



\k=l 
/ Nx 



{P,W) 



Y^i^lW^^M^) 



fc=l 
In (A (A)) 



(8b) 



The values of a and /(a) are obtained from the slopes of 
the linear fits of the averaged contributions in Eqs. ([5]) and 
([H|) (typical and ensemble respectively) versus In A, which 
is calculated for different values of the box-size I. Both 
averages in the thermodynamic limit are expected to tend 
towards the same spectrum for f{a) > 0. We refer the 
reader to Refs. [12, 13, 22] for a more detailed discussion 
of the relationship between both averages as used in an 
MFA. 



4 Partitioning into cubic boxes with integer 
ratios L/l 

The most simple way of partitioning the system with lin- 
ear size L is to use an isotropic cubic grid which can fit 
exactly the system. In this case the box-sizes in each di- 
rection satisfy Ix = ly = Iz with values of I such that 
L/l = n € N, and thus we always cover the system of 
size L^ with an integer number of boxes, as illustrated 
in Fig. Hal The singularity spectrum f{a) obtained using 




Fig. 2: Singularity spectrum obtained from box-size scal- 
ing of the typical average of Pq for L = 60 averaging over 
10"^ states using different schemes of box-partitioning: cu- 
bic boxes for integer L/l (solid black), cubic boxes for all 
values of / (grey), non-cubic boxes (dashed black). The 
symbols denote the position of g = 2 (left), q = (centre) 
and g = — 4 (right). 



this partitioning method for a system with i = 60 af- 
ter taking the typical average over 10'^ states is shown in 
Fig. [2] It must be emphasized that since we are dealing 
with a discrete system, the values of the box-size must be 
larger than the lattice spacing, in order to observe prop- 
erly the multifractal fiuctuations of the distribution [1,11]. 
Usually we consider values for the box-size in the interval 
I E [10, L/2]. The points used for the linear fits shown in 
Fig. El correspond to / = 10,12,15,20,30. Although easy 
to implement numerically, the drawback of the present 
method is that depending on the system size sometimes 
only a few values for I are allowed and that imposes a 
restriction on the reliability of the fits to obtain aq and 
fq. In the following sections we consider several addilSional 
partitioning strategies and compare their performance to 
the basic integer cubic-boxes technique. 



5 Using cubic boxes with unrestricted values 
for Z 

As illustrated in panel (b) of Fig. [U an alternative par- 
titioning method would be to use all values of I in the 
range 10 ^ / < L/2. The system is partitioned into cu- 
bic boxes without imposing the restriction of L/l being 
an integer. As such, there are values of I where some 
of the outer boxes will have a lack of sites. This can be 
overcome by using periodic boundary conditions where all 
necessary values \tpi\'^ will be repeated until all boxes are 
properly filled. One must realise that although the diago- 
nalisation process to obtain the eigenstates uses PBC to 
minimise edge effects, no true periodicity pattern exists 
in the wavefunction itself. However by using this parti- 
tioning scheme we are imposing a periodicity pattern on 
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Fig. 3: Linear fits of Eqs. © for a^^P (left) and L'yP (right) 
values of the singularity spectra shown in Fig.[2]using dif- 
ferent schemes of box-partitioning: cubic boxes for integer 
Ljl (black, square), cubic boxes for all values of / (grey, 
circle), non-cubic boxes (black, triangle). For the method 
of cubic boxes for all values of /, the data point correspond- 
ing to the average contribution for each A (plus symbol) 
is also shown. The dotted lines separate fits for q = 2 
(top) from g = -4 (bottom). The values of a'^P and /'^p 
are given by the slopes of the fits. Data points have been 
properly shifted vertically to ensure optimal visualisation. 
In all cases, the standard deviations are contained within 
symbol size. 



the distribution of iV'iP that in principle might distort its 
multifractal properties. Furthermore, we must emphasize 
here that repeating j^^p as seen in Fig. llbl will mean en- 
larging the system into L' and hence a renormalisation of 
the wavefunctions^ must be carried out if one wants to use 
the equations in Section [3l In Fig. [21 a comparison can 
be found between the typical singularity spectra obtained 
using the current box-partitioning method and the one in- 
volving cubic-boxes with integer Ljl. The corresponding 
linear fits used to obtain a^^P and /*^p for a couple of 
values of q are shown in Fig. [3l 

Using this partitioning approach the number of avail- 
able A values in the fits is not increased, since whenever 
the system is enlarged the correct value for the scaling 
variable must be also redefined to A = //£', which is the 
inverse of the number of boxes needed to cover the sys- 
tem in each dimesion. Therefore different box-sizes give 
contributions to the same A value in the fit. All these con- 
tributions as well as the corresponding weighted average 
are explicitly shown in Fig. [31 The linear fits performed 
to obtain the singularity spectrum are based on a x^ min- 
imisation taking into account the individual standard de- 
viations of each point. The first thing to notice is that the 
uncertainties of the individual contributions to a given A 
are clearly in conflict with the range of values spanned 



^ Alternatively the scaling law for the gIPR, as well as the 
definitions of T"(g), q, and f{ceq) can be generalised to account 
for unnormalised distributions. 



by the points. In fact we have observed that there exists 
a defined tendency of behaviour for these contributions. 
In the fits for aq (left panel in Fig. [3]), for large positive 
q and a fixed value of A the higher contribution is given 
by the box-size / which does not require an extension of 
the system. As I is increased and the system size is en- 
larged, preserving the value of A, the contributions are 
progressively smaller. This makes the value of the slope, 
that is aq, go systematically to smaller values at the end 
of the left branch of the spectrum, when compared to the 
method in Sec. [31 For high negative q the behaviour is 
the opposite, for a given A the additional contributions 
of the box-sizes requiring system enlargement grow with 
the box-size, and therefore the slope attains higher val- 
ues, i.e. aq moves towards higher values at the end of the 
right branch of the spectrum. This leads to a systematic 
broadening of /(a) as shown in Fig. [2l . The effect of the 
additional contributions in the fits for fq (right panel in 
Fig. [3]) causes a decrease of the slopes for high q, either 
negative or positive. We have also checked that the shape 
of the singularity spectrum is strongly dependent on the 
range of box-sizes taken into account for the linear fit. All 
these facts clearly make the present partitioning method 
very unreliable and not suitable to perform a numerical 
MFA. 



6 Partitioning with rectangular boxes 

Another strategy is to consider an anisotropic non-cubic 
box partitioning as shown in Fig. [Tel with at least one 
of the linear sizes lx,ly,h different from the other two. 
For simplicity we will consider values for the box sides 
such that L/lx.y^z is an integer. Hence in all directions the 
system can be covered exactly with an integer number of 
boxes, although this number can be different for each di- 
rection x,y,z. The scaling parameter used in the present 
case is A = l/L with an effective linear box-size defined as 
I ~ (Ixlylz)^^^- Let us note that different combinations of 
{Ixi ly,h} can lead to the same I and thus we could have 
different contributions to the same A, as in the method 
described in Sec.[5l In the fits of Fig. [3l only the averaged 
contribution for each A is shown. The corresponding sin- 
gularity spectrum (Fig. [2]) is similar to the /(a) obtained 
using only cubic boxes with integer L/l. We again note a 
tendency to broaden but it is less pronounced than when 
using unrestricted values for the box-size. Although in the 
present case the number of points used in the linear fits is 
noticeably increased as seen in Fig. [31 no significant im- 
provement with regards to the reliability of the f{a) has 
been observed. 



7 Using multiple origins for the box partitions 

An further possibility to increase the reliability of the mul- 
tifractal spectrum would be to consider different origins 
for the partitioning method as demonstrated in Fig. Ildl 
This means that instead of considering the box-partitioning 
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Fig. 4: Singularity spectrum obtained from box-size scal- 
ing for L = 60 averaging over 10^ states. Ensemble average 
using cubic boxes for integer L/l with single-origin (black) 
and multiple origins (grey) for the box partitioning. The 
error bars are equal to one standard deviation. The values 
for the linear correlation coefficient r^ and quality-of-fit 
parameter Q of the fits to obtain a™*^ and /™'' in the 
single origin and multiple origin cases are shown in the 
bottom shaded panels as labelled. 



from a single origin, we could use different points in the 
system as origins. This is equivalent, once we have done 
the partitioning for a given box-size I and we have a rigid 
partition grid, to shift the whole wavefunction at the same 
time one site at a time in all directions a maximum of I 
sites. Here, one uses the PBC in such a way that the wave- 
function sites that are left of the grid in one direction en- 
ter the system through the opposite end. Note that each 
value of iV'jP is used only once. Hence for each box-size 
we have l^ different non-equivalent origins for the box- 
partitioning. Therefore for a given box-size I the number 
of contributions in the average of the gIPR Pq is multiplied 
by l^ . Since the ensemble average is specially sensitive to 
the number of disorder realizations considered, we show 
in Fig. IHhow f™^{a) behaves when multiple origins are 
taken into account. The change with respect to the single- 
origin situation [12] is really insignificant, there is almost 
a perfect overlap for both lines. But the standard devi- 
ations are reduced almost to zero when multiple origins 
are considered. This is due to the fact that the standard 
deviations of the averaged contributions for the points in 
the fits, include a term 1/vW where M is the number 
of states we average over. Since every origin counts as a 
different contribution, then the uncertainties are greatly 
reduced. However, let us show that this reduction is in 
fact misleading and does not really mean an increase in 
reliability. In the bottom panels of Fig. 21 we compare the 
quality-of-fit parameter Q of the linear fits to obtain a and 
f{a) for the single-origin and multiple-origin techniques. 
In both cases the linear correlation coefficient r^ is close 
to 1 in the whole a range indicating a good linear be- 



haviour. But in the multiple-origin situation, the values of 
Qa and Qfia) decrease considerably suggesting that the 
uncertainties of the points in the fits have been clearly 
underestimated. Therefore there is no gain in reliability. 
Considering multiple origins maybe seen as considering 
eigenstates of other disordered systems which correspond 
to cyclic permutations in the three spatial directions of the 
initial disorder realization. Clearly, such transformations 
will not give independent disorder realizations. In fact the 
contribution of all these wavefunctions related by cyclic 
permutations must be very similar, since their probability 
distribution function for |-0|| values is identical. Therefore 
it is not correct to consider their contribution as a differ- 
ent disorder realization. That explains why the spectrum 
hardly changes when compared to the single-origin situa- 
tion. From our point of view this strategy should not be 
considered as a good method to reduce the uncertainty. It 
must also be emphasized that the multiple-origin strategy 
is very expensive in terms of computational time. 



8 An adaptive MFA-fit strategy 

Aside from the different schemes of partitioning, let us also 
study another strategy related to the selection of values of 
box-sizes that are considered for the linear fit. Having in 
mind that the proper region of values for I to do the scal- 
ing might also depend on the value of q, we implemented 
an adaptive-linear-fit strategy. For a given value of q, we 
consider a window with a certain number of points to do 
the fit and we maximise the values of r^ and Q by shifting 
this window throughout all the available interval of values 
for /, only for box-sizes with integer L/l. In this way for 
each q the best linear behaviour is achieved using a dif- 
ferent region of contiguous I values. In this case we have 
considered a system with size L ~ 210, for which 14 avail- 
able values of I exist in the interval I G [2,L/2], and we 
tried different widths of the window for /. Unfortunately 
using this strategy we have not been able to see proper 
multifractal spectra. Due to the fluctuations of the uncer- 
tainty for the points in the fits, sometimes the shifting of 
the window of /-values does not follow a smooth tendency 
with q, moreover it also happens that for a given q the 
best fit for aq and fq is achieved in a different range of I. 
These effects give rise to the appearance of discontinuities 
and irregularities in /(a). Apparently our initial premise 
is not true. 



9 Conclusion 

We have shown that the simple box-partitioning MFA 
method based on the use of cubic boxes with integer side 
length L/l and a single-origin for the partitioning, al- 
though bearing some limitations, is numerically the most 
reliable one. It is also optimal in terms of time and com- 
putational requirements, and it has already been success- 
fully implemented to obtain the scaling behaviour at the 
Anderson transition for very large system sizes [11,12]. 
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Therefore it should be considered the method of choice. 20. F. Milde, R. A. Romer, M. Schre iber, and V. Uski, E ur. 

Furthermore, let us point out that recently [23] an ana- Phys. J. B 15, 685 (2000), ArXiv: "cond-mat/9911029' 

lytical symmetry relationship for /(a) in the thermody- 21. M. BoUhofer and Y. Notay, Comp. Phys. Comm. 177, 951 

namic limit has been given. This relation states that the (2007). 

spectrum must be contained within ^ a ^ 6 for bulk 22. A. D. Mirlin and F. Evers, Phys. Rev. B 62, 7920 (2000). 

systems [24, 25]. Our resuhs in Fig. [2]show that the strat- 23. A. D. Mirlin, Y. V. Fyodorov, A. Mildenberg, and F. Evers, 

egy of cubic boxes with integer sizes L/l and the use of P^ys. Rev. Lett. 97, 046803 (2006). 

a single origin indeed gives the best agreement with the 24. H. Obuse et al, Phys. Rev. Lett. 98, 156802 (2007). 

syminetry [11, 12]. 25. H. Obuse et al, Physica E 40, 1404 (2008). 
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